* Webb, Linn an Lebo Procedure for independently-generated unit roots DGP

* Basics
cls 
cscript
set matsize 11000

* Set working directory
cd "/Users/alikagalwala/Library/CloudStorage/Dropbox/Whitten_Kagalwala/TS/JOP/Acceptance/KW_Replication/WLL Bounds Approach"

* Run the file with the simulation program
do "WLL-Coverage-DGP.do"

// Monte Carlo Analyses

* Create a matrix to store results
matrix mc_summary_stats = J(6, 23, .)

matrix colnames mc_summary_stats = t phi_y phi_x /// mc identifiers 
lrm_pa ftr_pa inc_pa rej_pa ///
lrm_gecm ftr_gecm inc_gecm rej_gecm ///
lrm_adl ftr_adl inc_adl rej_adl ///
lrm_adl22 ftr_adl22 inc_adl22 rej_adl22 ///
lrm_adl33 ftr_adl33 inc_adl33 rej_adl33

* Create a counter to store results in the matrix
local j = 1

foreach t in 25 50 75 150 500 1000 {
					parallel sim, expr(lrm_pa = r(lrm_pa) ftr_pa = r(ftr_pa) inc_pa = r(inc_pa) rej_pa = r(rej_pa) ///
					lrm_gecm = r(lrm_gecm) ftr_gecm = r(ftr_gecm) inc_gecm = r(inc_gecm) rej_gecm = r(rej_gecm) ///
					lrm_adl = r(lrm_adl) ftr_adl = r(ftr_adl) inc_adl = r(inc_adl) rej_adl = r(rej_adl) ///
					lrm_adl22 = r(lrm_adl22) ftr_adl22 = r(ftr_adl22) inc_adl22 =  r(inc_adl22) rej_adl22 = r(rej_adl22) ///
					lrm_adl33 = r(lrm_adl33) ftr_adl33 = r(ftr_adl33) inc_adl33= r(inc_adl33) rej_adl33 = r(rej_adl33)) ///
					reps(1000) noisily seeds(12344 12344 12344 12344 12344 12344 12344 12344): NG, obs(`t') phi_y(1) phi_x(1)
			
					matrix mc_summary_stats[`j', 1] = `t'
					matrix mc_summary_stats[`j', 2] = 1
					matrix mc_summary_stats[`j', 3] = 1

					* PA model
					qui su lrm_pa
					matrix mc_summary_stats[`j', 4] = r(mean)
					qui su ftr_pa
					matrix mc_summary_stats[`j', 5] = r(mean)
					qui su inc_pa
					matrix mc_summary_stats[`j', 6] = r(mean)
					qui su rej_pa
					matrix mc_summary_stats[`j', 7] = r(mean)
					
					* GECM model
					qui su lrm_gecm
					matrix mc_summary_stats[`j', 8] = r(mean)
					qui su ftr_gecm
					matrix mc_summary_stats[`j', 9] = r(mean)
					qui su inc_gecm
					matrix mc_summary_stats[`j', 10] = r(mean)
					qui su rej_gecm
					matrix mc_summary_stats[`j', 11] = r(mean)
			
					* ADL model
					qui su lrm_adl
					matrix mc_summary_stats[`j', 12] = r(mean)
					qui su ftr_adl
					matrix mc_summary_stats[`j', 13] = r(mean)
					qui su inc_adl
					matrix mc_summary_stats[`j', 14] = r(mean)
					qui su rej_adl
					matrix mc_summary_stats[`j', 15] = r(mean)
					
					* ADL(2,2) Model
					qui su lrm_adl22
					matrix mc_summary_stats[`j', 16] = r(mean)
					qui su ftr_adl22
					matrix mc_summary_stats[`j', 17] = r(mean)
					qui su inc_adl22
					matrix mc_summary_stats[`j', 18] = r(mean)
					qui su rej_adl22
					matrix mc_summary_stats[`j', 19] = r(mean)
					
					* ADL(3,3) Model
					qui su lrm_adl33
					matrix mc_summary_stats[`j', 20] = r(mean)
					qui su ftr_adl33
					matrix mc_summary_stats[`j', 21] = r(mean)
					qui su inc_adl33
					matrix mc_summary_stats[`j', 22] = r(mean)
					qui su rej_adl33
					matrix mc_summary_stats[`j', 23] = r(mean)
					* Iterate the counter
					local j = `j' + 1
		}

// Store MC matrix results as a dataframe
clear
set obs 6
svmat double mc_summary_stats, names(col)
save bounds_coverage.dta, replace
